C
C     ..................................................................
C
C        SUBROUTINE GELB
C
C        PURPOSE
C           TO SOLVE A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH A
C           COEFFICIENT MATRIX OF BAND STRUCTURE.
C
C        USAGE
C           CALL GELB(R,A,M,N,MUD,MLD,EPS,IER)
C
C        DESCRIPTION OF PARAMETERS
C           R      - M BY N RIGHT HAND SIDE MATRIX (DESTROYED).
C                    ON RETURN R CONTAINS THE SOLUTION OF THE EQUATIONS.
C           A      - M BY M COEFFICIENT MATRIX WITH BAND STRUCTURE
C                    (DESTROYED).
C           M      - THE NUMBER OF EQUATIONS IN THE SYSTEM.
C           N      - THE NUMBER OF RIGHT HAND SIDE VECTORS.
C           MUD    - THE NUMBER OF UPPER CODIAGONALS (THAT MEANS
C                    CODIAGONALS ABOVE MAIN DIAGONAL).
C           MLD    - THE NUMBER OF LOWER CODIAGONALS (THAT MEANS
C                    CODIAGONALS BELOW MAIN DIAGONAL).
C           EPS    - AN INPUT CONSTANT WHICH IS USED AS RELATIVE
C                    TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.
C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
C                    IER=0  - NO ERROR,
C                    IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAME-
C                             TERS M,MUD,MLD OR BECAUSE OF PIVOT ELEMENT
C                             AT ANY ELIMINATION STEP EQUAL TO 0,
C                    IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
C                             CANCE INDICATED AT ELIMINATION STEP K+1,
C                             WHERE PIVOT ELEMENT WAS LESS THAN OR
C                             EQUAL TO THE INTERNAL TOLERANCE EPS TIMES
C                             ABSOLUTELY GREATEST ELEMENT OF MATRIX A.
C
C        REMARKS
C           BAND MATRIX A IS ASSUMED TO BE STORED ROWWISE IN THE FIRST
C           ME SUCCESSIVE STORAGE LOCATIONS OF TOTALLY NEEDED MA
C           STORAGE LOCATIONS, WHERE
C             MA=M*MC-ML*(ML+1)/2    AND    ME=MA-MU*(MU+1)/2    WITH
C             MC=MIN(M,1+MUD+MLD),  ML=MC-1-MLD,  MU=MC-1-MUD.
C           RIGHT HAND SIDE MATRIX R IS ASSUMED TO BE STORED COLUMNWISE
C           IN N*M SUCCESSIVE STORAGE LOCATIONS. ON RETURN SOLUTION
C           MATRIX R IS STORED COLUMNWISE TOO.
C           INPUT PARAMETERS M, MUD, MLD SHOULD SATISFY THE FOLLOWING
C           RESTRICTIONS     MUD NOT LESS THAN ZERO
C                            MLD NOT LESS THAN ZERO
C                            MUD+MLD NOT GREATER THAN 2*M-2.
C           NO ACTION BESIDES ERROR MESSAGE IER=-1 TAKES PLACE IF THESE
C           RESTRICTIONS ARE NOT SATISFIED.
C           THE PROCEDURE GIVES RESULTS IF THE RESTRICTIONS ON INPUT
C           PARAMETERS ARE SATISFIED AND IF PIVOT ELEMENTS AT ALL
C           ELIMINATION STEPS ARE DIFFERENT FROM 0. HOWEVER WARNING
C           IER=K - IF GIVEN - INDICATES POSSIBLE LOSS OF SIGNIFICANCE.
C           IN CASE OF A WELL SCALED MATRIX A AND APPROPRIATE TOLERANCE
C           EPS, IER=K MAY BE INTERPRETED THAT MATRIX A HAS THE RANK K.
C           NO WARNING IS GIVEN IF MATRIX A HAS NO LOWER CODIAGONAL.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           SOLUTION IS DONE BY MEANS OF GAUSS ELIMINATION WITH
C           COLUMN PIVOTING ONLY, IN ORDER TO PRESERVE BAND STRUCTURE
C           IN REMAINING COEFFICIENT MATRICES.
C
C     ..................................................................
C
      SUBROUTINE GELB(R,A,M,N,MUD,MLD,EPS,IER)
C
C
      DIMENSION R(1),A(1)
C
C     TEST ON WRONG INPUT PARAMETERS
      IF(MLD)47,1,1
    1 IF(MUD)47,2,2
    2 MC=1+MLD+MUD
      IF(MC+1-M-M)3,3,47
C
C     PREPARE INTEGER PARAMETERS
C        MC=NUMBER OF COLUMNS IN MATRIX A
C        MU=NUMBER OF ZEROS TO BE INSERTED IN FIRST ROW OF MATRIX A
C        ML=NUMBER OF MISSING ELEMENTS IN LAST ROW OF MATRIX A
C        MR=INDEX OF LAST ROW IN MATRIX A WITH MC ELEMENTS
C        MZ=TOTAL NUMBER OF ZEROS TO BE INSERTED IN MATRIX A
C        MA=TOTAL NUMBER OF STORAGE LOCATIONS NECESSARY FOR MATRIX A
C        NM=NUMBER OF ELEMENTS IN MATRIX R
    3 IF(MC-M)5,5,4
    4 MC=M
    5 MU=MC-MUD-1
      ML=MC-MLD-1
      MR=M-ML
      MZ=(MU*(MU+1))/2
      MA=M*MC-(ML*(ML+1))/2
      NM=N*M
C
C     MOVE ELEMENTS BACKWARD AND SEARCH FOR ABSOLUTELY GREATEST ELEMENT
C     (NOT NECESSARY IN CASE OF A MATRIX WITHOUT LOWER CODIAGONALS)
      IER=0
      PIV=0.
      IF(MLD)14,14,6
    6 JJ=MA
      J=MA-MZ
      KST=J
      DO 9 K=1,KST
      TB=A(J)
      A(JJ)=TB
      TB=ABS(TB)
      IF(TB-PIV)8,8,7
    7 PIV=TB
    8 J=J-1
    9 JJ=JJ-1
C
C     INSERT ZEROS IN FIRST MU ROWS (NOT NECESSARY IN CASE MZ=0)
      IF(MZ)14,14,10
   10 JJ=1
      J=1+MZ
      IC=1+MUD
      DO 13 I=1,MU
      DO 12 K=1,MC
      A(JJ)=0.
      IF(K-IC)11,11,12
   11 A(JJ)=A(J)
      J=J+1
   12 JJ=JJ+1
   13 IC=IC+1
C
C     GENERATE TEST VALUE FOR SINGULARITY
   14 TOL=EPS*PIV
C
C
C     START DECOMPOSITION LOOP
      KST=1
      IDST=MC
      IC=MC-1
      DO 38 K=1,M
      IF(K-MR-1)16,16,15
   15 IDST=IDST-1
   16 ID=IDST
      ILR=K+MLD
      IF(ILR-M)18,18,17
   17 ILR=M
   18 II=KST
C
C     PIVOT SEARCH IN FIRST COLUMN (ROW INDEXES FROM I=K UP TO I=ILR)
      PIV=0.
      DO 22 I=K,ILR
      TB=ABS(A(II))
      IF(TB-PIV)20,20,19
   19 PIV=TB
      J=I
      JJ=II
   20 IF(I-MR)22,22,21
   21 ID=ID-1
   22 II=II+ID
C
C     TEST ON SINGULARITY
      IF(PIV)47,47,23
   23 IF(IER)26,24,26
   24 IF(PIV-TOL)25,25,26
   25 IER=K-1
   26 PIV=1./A(JJ)
C
C     PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R
      ID=J-K
      DO 27 I=K,NM,M
      II=I+ID
      TB=PIV*R(II)
      R(II)=R(I)
   27 R(I)=TB
C
C     PIVOT ROW REDUCTION AND ROW INTERCHANGE IN COEFFICIENT MATRIX A
      II=KST
      J=JJ+IC
      DO 28 I=JJ,J
      TB=PIV*A(I)
      A(I)=A(II)
      A(II)=TB
   28 II=II+1
C
C     ELEMENT REDUCTION
      IF(K-ILR)29,34,34
   29 ID=KST
      II=K+1
      MU=KST+1
      MZ=KST+IC
      DO 33 I=II,ILR
C
C     IN MATRIX A
      ID=ID+MC
      JJ=I-MR-1
      IF(JJ)31,31,30
   30 ID=ID-JJ
   31 PIV=-A(ID)
      J=ID+1
      DO 32 JJ=MU,MZ
      A(J-1)=A(J)+PIV*A(JJ)
   32 J=J+1
      A(J-1)=0.
C
C     IN MATRIX R
      J=K
      DO 33 JJ=I,NM,M
      R(JJ)=R(JJ)+PIV*R(J)
   33 J=J+M
   34 KST=KST+MC
      IF(ILR-MR)36,35,35
   35 IC=IC-1
   36 ID=K-MR
      IF(ID)38,38,37
   37 KST=KST-ID
   38 CONTINUE
C     END OF DECOMPOSITION LOOP
C
C
C     BACK SUBSTITUTION
      IF(MC-1)46,46,39
   39 IC=2
      KST=MA+ML-MC+2
      II=M
      DO 45 I=2,M
      KST=KST-MC
      II=II-1
      J=II-MR
      IF(J)41,41,40
   40 KST=KST+J
   41 DO 43 J=II,NM,M
      TB=R(J)
      MZ=KST+IC-2
      ID=J
      DO 42 JJ=KST,MZ
      ID=ID+1
   42 TB=TB-A(JJ)*R(ID)
   43 R(J)=TB
      IF(IC-MC)44,45,45
   44 IC=IC+1
   45 CONTINUE
   46 RETURN
C
C
C     ERROR RETURN
   47 IER=-1
      RETURN
      END
